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We consider open quantum systems with dynamics described by master equations that have 
perturbative expansions in the system-environment interaction. We show that, contrary to intuition, 
full-time solutions of order-2n accuracy require an order-(2n + 2) master equation. We give two 
examples of such inaccuracies in the solutions to an order-2n master equation: order-2n inaccuracies 
in the steady state of the system and order-2n positivity violations, and we show how these arise 
in a specific example for which exact solutions are available. This result has a wide-ranging impact 
on the validity of coupling (or friction) sensitive results derived from second-order convolutionless, 
Nakajima-Zwanzig, Redfield, and Born-Markov master equations. 



I. INTRODUCTION 

An open quantum system is a quantum system that in- 
teracts with some environment whose degrees of freedom 
have been coarse grained over (i.e., traced out), and its 
dynamics are described by a master equation governing 
the reduced density matrix p. Exact master equations for 
the stochastic dynamics of open quantum systems are, in 
general, out of reach. However, arbitrary-order perturba- 
tive master equations (in the system-environment inter- 
action) can be derived in a variety of different ways [1-3] 
and find application in many branches of physics and 
chemistry [4-7]. In the time-local representation (also 
called the convolutionless or Markovian representation), 
the dynamics of the reduced density matrix of the system 
can be expressed with a quantum Liouville equation 

^p(t)=C(t)p(t), (1) 

where, despite the apparent time-local form, non- 
Markovian behavior may be encapsulated in the time 
dependence of the Liouvillian C(t). As a perturbative 
approximation, C(t) is expanded in powers of the system- 
environment interaction, c, and truncated to some order. 

We will consider perturbative master equations where 
the Liouvillian C{t) is time independent at zeroth order 
and asymptotically constant for late times. We will as- 
sume that the perturbative expansion of C(t) is in even 
powers of the coupling, because as we explain in the 
next section, this can naturally arise from a microscopic 
derivation of the open-system dynamics. The expansion 
of C(t) will then take the form 

oo 

C(t) = J2^[2n](t), (2) 
n=0 

C m {t)p = [-iU,p], (3) 

where £[ 2 „](£) = 0(c 2n ) and to zeroth order the system 
is driven in a unitary manner by its Hamiltonian H. 

The most well-known perturbative master equation is 
the second-order master equation, as it can be equiva- 
lent to the Redfield and Born-Markov master equations. 



This is partly due to the fact that in the Markovian limit, 
the second-order master equation is exact. But equiva- 
lence with the previous approximate master equations 
does not carry to fourth order and there perturbation 
theory is strictly superior. One might easily assume that 
solving the second-order master equation defined by the 
Liouvillian £[ ] + Cm\ would yield a solution that would 
match the exact solution to the exact master equation 
up to second order, having error terms of order C(c 4 ); 
however we will show that in general they will differ by 
second-order terms, so that one can only say they are in 
perturbative agreement at zeroth order. 

One very significant implication of these facts is for 
positivity. Not being exact, nor generally of Lindblad 
form [8, 9], a perturbative master equation is not guar- 
anteed to yield a dynamical map with exact complete 
positivity. Solutions can and should be completely pos- 
itive to the relevant perturbative order, and as we show 
in this work that order is not what one might naively ex- 
pect. Solutions to the second-order master equation can 
violate positivity by an amount that is C(c 2 ). We show 
that to find solutions good to second-order, canonical 
perturbation theory generally demands the fourth-order 
Liouvillian. 

In our penultimate section, we use the example of 
quantum Brownian motion to show how naive usage of 
the second-order master equation leads to an underes- 
timation of the second-order position uncertainty of a 
damped quantum oscillator (such as a nanomechanical 
resonator). We would refer the reader to this section to 
see an immediate physical implication of the accuracy 
loss we describe in this work. 

Finally, in our core proof we employ the formalism of 
canonical perturbation theory as applied to master equa- 
tions [10]. If the reader is interested in seeing a specific 
example problem worked out in its entirety with these 
techniques, we would refer the reader to Ref. [ ] which 
considered neutral atoms held fixed in the electromag- 
netic field. There it was found that the second-order mas- 
ter equation is completely inadequate in assessing sudden 
death of entanglement at low temperature. It was also in 
this work that we noticed the relation between second- 
order positivity loss and low temperature, which will be 
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discussed later. 



II. MASTER EQUATIONS FROM A 
MICROSCOPIC MODEL 

In a system derived from a microscopic model, the 
coarse-grained environment can act as a source of noise, 
dissipation, and decoherence; thus its influence provides 
a model of dissipative quantum mechanics more general 
than Markovian (white-noise) models, which can be con- 
structed phenomenologically. Such microscopic models 
can still lead to a stationary Liouvillian in the late-time 
limit, even when the dynamics they describe are non- 
Markovian. 

Given a stationary system Hamiltonian and station- 
ary bath correlations, Gaussian noise distributionals (e.g. 
noise generated via linear coupling to an environment of 
harmonic oscillators) may allow the master equation to 
have a stationary late-time limit [10] 

C(oo) = lim C(t) , (4) 

t— >oo 

so that the late-time and weak-coupling limits commute; 
otherwise perturbation theory cannot be used for long 
durations of time. Gaussian noise processes are catego- 
rized by their second-order noise correlation, and whether 
or not the master equation will have a stationary limit 
is dependent upon how localized this noise correlation is. 
Well-localized noise correlations (e.g. Gaussian or expo- 
nential) can lead to a very well-behaved master equation, 
whereas long-ranged noise correlations (e.g. Cauchy) can 
produce a more pathological master equation which can- 
not be accurately analyzed in a perturbative fashion. Ex- 
act examples of this phenomena are given in Ref. [12] in 
the context of quantum Brownian motion with Ohmic 
and sub-Ohmic couplings. Moreover, the exact solutions 
p(t) can be very well-behaved even if C(t) is not. Marko- 
vian representations (and, more generally, effective equa- 
tions of motion) are not always suitable. 

III. INDETERMINACY OF SOLUTIONS 

Before determining what the appropriate level of ac- 
curacy is for the solutions of perturbative master equa- 
tions, we will first demonstrate that there is issue with 
the naive expectation of order-2n accuracy. This argu- 
ment is a generalization of one found in Ref. [13], where 
the discrepancy was noticed for the second-order equilib- 
rium state. Let P(2„)(i) be any solution which satisfies 
the master equation (and is supposedly accurate) to or- 
der 2n, then 

j t P {2n) (t)=C(t)p {2n) (t) + 0(c^). (5) 
Furthermore consider the order-2n state 

P'(2n) (t) = P(2n) (*) + S P[2n] W > (6) 



where 6p^ 2n ] is an order-2n traceless and diagonal (in the 
energy basis) perturbation for which 

6 P[2n] (0) = 0, (7) 

j t 6p l2n] (t) = 0(c 2n+2 ), (8) 

so that both P( 2n ) W an d P(2n) W share the same ini- 
tial conditions, and the discrepancy <5p[2n] (t) grows slowly 
with the perturbation as to also satisfy 

j t p{ 2n) (t)=C(t)p' {2n) (t)+0(c^). (9) 

given that Co <5pr 2n i (t) = by construction. This demon- 
strates that, for non-perturbative durations of time, there 
is an order 2n ambiguity in the stationary (e.g. diagonal) 
entries of all solutions if one only compares terms up to 
order 2n. This proof also applies to time-nonlocal mas- 
ter equations, replacing perturbative contributions to the 
Liouvillian with corresponding memory-kernel operators. 
Next we will proceed to our main proofs where we show 
how this issue arises, that this is the full extent of the 
problem, and precisely how it can be remedied. 

IV. LATE-TIME ACCURACY 

It is clear that if Eqs. (1) and (2) are well defined then 
for sufficiently short times an order-2n master equation 
(in which the sum in Eq. (2) only includes terms up to or- 
der 2n) can produce a solution that is also accurate to or- 
der 2n. We find that for longer spans of time, and in par- 
ticular the late-time regime wherein the master equation 
assumes its stationary limit (4), solutions to the order-2n 
master equation are only accurate to order 2n — 2. The 
reason is an ultimately mundane but slightly subtle result 
of degenerate perturbation theory. In this section we will 
address the late-time stationary dynamics, and then in 
following sections we will address the full-time dynamics, 
including the crossover from consistent accuracy to loss 
of accuracy. 

Assuming we have the perturbative expansion of a sta- 
tionary master equation (i.e., an expansion of £), we then 
seek perturbative solutions obtained by applying canon- 
ical perturbation theory of the eigenvalue problem [10] 

£o = /o, (10) 

where o is a Hilbert-space eigen-operator and / is its 
corresponding eigen-value. We already know the zeroth- 
ordcr solutions 

C[ ] \Ui)(uJj\ = -IWij \Ui)(u]j\ , (11) 

where H \u>) = uj \uj) and denote the (free) 

energy basis of the system. In the appropriate regime of 
validity, exact solutions to the perturbative master equa- 
tion should agree with the perturbative solutions to the 
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exact master equation up to the appropriate order. Note 
that perturbation theory with master equations is always 
degenerate perturbation theory as uia = LUjj = 0. This 
inevitably-degenerate subspace corresponds to the space 
of operators that are diagonal in the energy basis of the 
free system. For simplicity let us assume no other degen- 
eracy in the spectrum of the free Liouvillian (though the 
possibility of extra degeneracy or near degeneracy arising 
from resonance can be suitably dealt with). 

Perturbation theory tells us that the second-order cor- 
rections to all eigenvalues and eigenoperators of C outside 
the degenerate subspace (off-diagonal operators) can be 
computed using only the second-order master equation: 



J ij 



{u)e\o§ \w r ) 



(Wi|.C[2]{|cJiXWj|} 
(u>i>\C[2]{\Ui}{Uj\}\Wj' 



-z(Wj 



) 



(12) 
(13) 



As is usual in degenerate perturbation theory, to com- 
pute corrections to eigen-operators from the degenerate 
subspace, which all satisfy Cw\ o' ' = 0, we must diago- 
nalize C in the degenerate subspace. This is equivalent to 
finding the correct linear combination of eigen-operators 
which branch under perturbation. The associated char- 
acteristic equation can be written 



Wo = /o, 

{6\i = (u>i\o\u>i) , 



(14) 
(15) 



where o denotes the degenerate-subspace projection of o 
represented as a vector, i.e. diagonal entries of the eigen- 
operator while in the free energy basis, and W defines 
the Pauli master equation 



(16) 



which is the degenerate-subspace projection of C repre- 
sented as a matrix, i.e. master-equation super-operators 
which map diagonal entries to diagonal entries. There- 
fore Eq. (14) must be solved for with Wpi exactly, and 
then the further effects of W[ 4 ] , W[ 6 ] , etc. , can be incorpo- 
rated via canonical perturbation theory. [Note that this 
is slightly more complicated than the usual canonical per- 
turbation in the Schrodinger equation where one knows 
the Hamiltonian perturbation exactly] The eigenval- 
ues obtained in diagonalizing Wr 2 i give the second-order 
corrections to the eigenvalues of C and the correct 
zcroth-order eigenoperators for the degenerate sub- 
space. Degenerate perturbation theory tells us that, in 
order to calculate each for the degenerate subspace, 
one actually requires Wmi from the fourth-order master 
equation; it will contribute the second-order correction 



Wr 



[4] 



[2] 



(17) 



where o* is the left eigen- vector of W such that o* W 



o* fi and o* 6 l 



Such corrections would be fourth 



order in a non-degenerate problem, but because the free 
Liouvillian is always degenerate, they become second or- 
der as the relevant lowest-order nonvanishing eigenvalue 
splitting is always second order here. Without this infor- 
mation from the fourth-order master equation, one can- 
not generate the complete second-order solution. 

Finally note that this requirement must extend even to 
exact solutions of the perturbative master equation. A 
perturbative solution to the second-order master equa- 
tion will be equivalent to solving the full master equa- 
tion perturbatively and then artificially setting £r 4 i and 
all higher-order contributions to the Liouvillian to van- 
ish. From this and the preceding perturbative analysis 
we know that the second-order perturbative solutions to 
the exact and second-order master equations must dif- 
fer by a term that is 0(c 2 ). Since the exact solutions 
to each given master equation must differ from the cor- 
responding second-order perturbative solutions by terms 
of C(c 4 ), we can conclude from our analysis that even the 
exact solution to the second-order master equation differs 
from the exact solution to the full master equation by a 
term of 0(c 2 ). In the final section we use the example 
of quantum Brownian motion, where an exact solution is 
available, to show that the second-order corrections aris- 
ing from the fourth-order Liouvillian are indeed present. 

More generally, the same argument tells us that while 
the short-time accuracy of an order-2n master equation 
can also be order 2n, the long-time accuracy can only be 
order 2n — 2. To obtain order-2n solutions one requires 
not only the order-2n master equation but in addition the 
order-(2n+2) Pauli master equation. In particular, the 
second-order master equation after taking the rotating- 
wave approximation [l ] will contain just enough terms 
to generate solutions which are accurate to zeroth-order 
[ !:]. The full second-order master equation improves 
upon this but not enough to generate the full second- 
order solutions. 

Among the information missing due to the second- 
order errors of the solution to the second-order master 
equation are important contributions to the asymptotic 
state of the system. When coupled to a thermal reservoir 
the system must asymptote to p oc e~^ H for vanishing 
system-environment coupling (though this may also hap- 
pen in other, very specific approximations [ ]). One of- 
ten desires to find the additional environmentally induced 
system-system correlations (and possibly entanglement) 
provided by perturbative corrections, but these will not 
be given correctly by directly finding the steady state 
of the second-order master equation. However, at least 
for zero-temperature noise, it is still possible to easily 
construct via other methods the order-2ro corrections us- 
ing only order-2n master equation coefficients and limits 
thereof [10, 11]. 

Another important characteristic that is mangled by 
the second-order master equation is positivity, as was 
mentioned in the introduction. The second-order inac- 
curacies that arise from using the second-order master 
equation imply that the diagonal elements of the density 
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matrix in the (free) energy basis are off by second-order 
terms. This can lead to second-order violations of posi- 
tivity. In fact, this is almost guaranteed at low tempera- 
ture, where any off-diagonal perturbation to the ground 
state will immediately cause second-order positivity vio- 
lation, given that the necessary inequality 

Pa Pjj > Pij Pji , (18) 

cannot be satisfied with the left-hand side vanishing at 
zcroth-order and not perturbed to the correct second- 
order values. 

V. FULL-TIME ACCURACY 

In analyzing the full-time accuracy of time-dependent 
master equations, first we will show that the short-time 
solutions are accurate to the order of the master equa- 
tion, and then we will show that longer-time solutions 
display accuracy loss. The timescale for this transition is 
determined by the frequency perturbations, e.g. l//[2l- 

To analyze the short-time behavior we rotate to the 
interaction picture defined 

p{t) = G^\t)p{t), (19) 
G Q (t)p = e- lHt pe +lHt , (20) 

wherein the master equation is now given by 

j t p(t)=5C(t)p(t), (21) 

dC(t) = G^(t)5C(t)G (t), (22) 
SC(t) = C{t) - C , (23) 

and so the interaction-picture dynamics are strictly per- 
turbativc. Short-time solutions can be obtained from the 
Neumann series 

p(t) = G(t)p(0), (24) 

G(t) = 1 + f dr6C(T) + [ dr [ dr' 8C{t) 8C{t') + ■■■ , 
Jo Jo Jo 

(25) 

where the order-2n solution is fully determined by 
£[2n](t)- However, such solutions are inherently secular 
in time. If /p] denotes the second-order frequency per- 
turbations, e.g. dissipation and diffusion rates, then the 
above solutions (at second order) are only good for times 
t <C 1/ /[2] • This is the regime wherein perturbative mas- 
ter equations are ensured to provide matching accuracy 
in their solutions. 

For longer spans of time, one must resort to time- 
ordered integration for solutions. For weak coupling the 
master equation can asymptote to its stationary value 
within timescales much shorter than l//[2], and so one 
can apply the stationary master equation and our corre- 
sponding proof of accuracy loss. More generally one may 



consider the behavior of the time-dependent eigen-value 
equation 

C(t)o(t) = f(t)o(t), (26) 

so that the time-translation generator may be given by 
its spectral decomposition 

C(t) = J2fk(t)o k (t)ot(t). (27) 

k 

Again, the order-2n master equation can only deter- 
mine the perturbatively-stationary eigen-operators o(t) 
to within order In — 2. Given that the time-dependent 
basis of the time-translation generator cannot be deter- 
mined to second order, neither can the solutions. 

One might be concerned with how the proof of short- 
time accuracy is compatible with this proof of full-time 
accuracy loss. In fact, the short-accuracy occurs within 
a span of time < t <C l//[2], which is not sufficient 
enough to accumulate full-order contributions from the 
perturbation. Therefore the regime of short-time accu- 
racy is a rather trivial result. 

VI. TIME NON-LOCAL ACCURACY 

Corresponding to the time-local master equation (1) is 
the time-nonlocal master equation 

j f p{t) = JdTK(t-T)p(r), (28) 

first derived via the projection-operator formalism of 
Nakajima [16] and Zwanzig [17]. The two representations 
are contrasted in Ref. [10, 18, 19]. The nonlocal kernel 
7C(i) also has a perturbative expansion with zeroth-order 
dynamics given by 

K [0] (t)=2S(t)C [0] . (29) 

which is time-local and unitary. Solutions are most easily 
calculated in the Laplace domain wherein one has the 
kernel 

/>oo 

K(s)= / dte- ts K.{t), (30) 
Jo 

£[0](s) =£[o] • (31) 

Perturbative solutions can then be acquired by solving 
the nonlocal eigen-value equation [10] 

it(s) 6(s) = k(s) 6(s) , (32) 

where from Eq. (31) the nonlocal eigen-system must be 
a perturbation of the free system-energy eigen-system, 
and therefore our proof of accuracy loss will carry over. 
The order-2n master equation can only determine the 
perturbatively-stationary eigen-operators 6(s) to within 
order 2n— 2. 
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VII. EXAMPLE: QBM 

As an example of an exactly-solvable open system, let 
us consider the master equation of an oscillator bilinearly 



coupled (position-position) to an environment of oscilla- 
tors initially in a thermal state [20]: 



dt 



p = [-t H R , p] - iT [x, {p, p}} - M D pp [x, [x, p}] - D xp [x, [p, p] 



(33) 



r 



where Hr is the system Hamiltonian but with frequency 
Or , r is the dissipation coefficient, D pp and D xp are the 
regular and anomalous diffusion coefficients. This mas- 
ter equation describes the dynamics of damped nano- 
mechanical resonators at low temperature, among other 
physical systems. 

In Ref. [12] exact solutions are given with full time de- 
pendence, and it is from this reference that we take all 
of the following results. Let us consider Ohmic coupling 
to the bath with damping kernel ^(t) = 2jq 5\(t), where 
5 a (t) is a representation of the delta function in the high- 
frequency cutoff limit A — > oo. [The damping kernel, and 
thus 70, is second order in the system-environment inter- 
action c] The homogeneous coefficients quickly asymp- 
tote to Or = O and r = 7q within the cutoff timescale, 
whereas the diffusion coefficients asymptote to 



D xp = +7 Im[2o] , 



2-f T + 7 Im I 




(34) 
(35) 



(37) 



mostly within the system timescale, but also hastened 
by temperature. In all coefficients we have neglected 
terms of order C(I/A). H here is the harmonic num- 
ber function, which is asymptotically logarithmic and yet 
H(0) = 0. Therefore both diffusion coefficients contain 
logarithmic cutoff sensitivities, though the sensitivity is 
present in the anomalous diffusion coefficient at second 
order, whereas it does not appear in the regular diffusion 
coefficient until fourth order. 

In the stationary limit, the system relaxes into a Gaus- 
sian state with phase-space covariancc 



v 2T PP 





D xp ) 



M n 

2T LJ PP 



(38) 



One can see that for a second-order master equation, the 
contribution from the regular diffusion D pp /T starts at 
zcroth order, while the contribution from anomalous dif- 
fusion D xp starts at second order. The full second-order 
contribution from the regular diffusion actually requires 
the fourth-order coefficients. 



In the exact calculation, or in any consistent per- 
turbative calculation, the logarithmic cutoff sensitivities 
present in the diffusion coefficients actually cancel in the 
position uncertainty. In this sense the anomalous diffu- 
sion coefficient acts as an anti-diffusion coefficient and 
this behavior will also occur for supra-Ohmic couplings. 
If one were to naively apply the second-order diffusion 
coefficients, and solve the master equation exactly, then 
one would obtain a mixed-order result and the logarith- 
mic cutoff sensitivities would not precisely cancel. The 
position uncertainty would contain a second-order neg- 
ative log (A) contribution. For sufficiently large cutoff 
frequencies, the Heisenberg uncertainty principle would 
be violated. For even larger frequencies, the covariancc 
would become negative. In any case the second-order 
master equation would produce a (supposedly) second- 
order position uncertainty which is an underestimation 
of the true second-order uncertainty. 



VIII. DISCUSSION 

We have shown that even when provided with a sta- 
tionary master equation describing dynamics that are 
amenable to perturbative solution, the solutions to an 
order-2n perturbative master equation are, in general, 
only accurate to order-(2n — 2), a step down from that 
of the master equation itself. This has a wide-range of 
implications upon the common use of second-order mas- 
ter equations and related master equations derived from 
second-order dynamics: the Redfield, Born-Markov, and 
many Lindblad equations. Moreover, not even a nonlo- 
cal representation, such as with the Nakajima-Zwanzig 
master equation can avoid this effect. This is to be ex- 
pected as a thorough analysis of time-local and nonlocal 
dynamics shows their asymptotics to be perturbatively 
the same [10]. 

To be more specific, the second-order master equation 
can provide all second-order timescales and off-diagonal 
density matrix elements (in the free energy basis). How- 
ever it can only provide the diagonal matrix elements 
with zeroth-order accuracy, and the missing information 
is the most relevant to positivity in the low-temperature 
regime. Therefore the second-order master equation can 
produce second-order positivity violations, whereas the 
full second-order solutions are positive to second-order. 
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Likewise, the steady state of the second-order master 
equation may only agree with the steady state of the 
full master equation to zeroth order. More generally, the 
predicted expectation of observables will typically be off 
by a second-order amount, and this certainly includes the 
energy or other quantities that were conserved at zeroth 
order. We have shown that this inaccuracy manifests it- 
self in the case of quantum Brownian motion through an 
underestimation of the position undertainty stationary 
limit. The same issue also affects the predicted dynam- 
ics for a collection of atoms interacting with a shared 
field [11], where the complete second-order solution is re- 
quired make correct predictions about the sudden death 
of entanglement. 

There are three mathematical limits in which the 
second-order master equation will give solutions accu- 
rate to second order: The first is early times, where t 
is small compared to any of the second-order damping 



time scales. The second is the Markovian limit, because 
in this limit the second-order master equation is exact. 
The third is the limit employed by Davies [21] where one 
rewrites the master equation in terms of the rescaled time 
parameter r = c 2 i and then takes the limit c — > (for 
t ^ this effectively amounts to taking a simultaneous 
t — > oo limit). In this limit all corrections to the eigen- 
operators of the Liouvillian vanish, and the only effect of 
the environment is to introduce damping rates through 
corrections to the eigenvalues, which are correctly cap- 
tured by the perturbative master equation. Thus, the 
inaccuracies of second-order master equation we have ad- 
dressed may be sufficiently suppressed even at late times 
if a physical system is sufficiently close to being described 
by one of these limits. Therefore our results should be 
most heeded in the non-Markovian regime of low tem- 
perature or long-ranged correlations in contexts where 
C(c 2 ) discrepancies are not negligible. 
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